Pregunta 1

2. 50 observaciones de una \(\mathcal{N}(0,1)\) y otras 50 observaciones de una \(\mathcal{N}(2,1)\). ¿Cómo se verán las 100 caras de Chernoff si \(X\) y \(Y\) definen la línea de la cara y la oscuridad del cabello?, ¿esperan caras similares?, ¿cuántas caras lucen como oservaciones de \(Y\) cuando aún son observaciones de \(X\)?

  ind      = matrix(0, ncol = 36)   # define una indicadora para el argumento which.row
  ind[,13] = 1                      # linea derecha de la cara
  ind[,14] = 1                      # oscuridad del cabello lado derecho
  ind[,31] = 1                      # linea izquierda de la cara
  ind[,32] = 1                      # oscuridad izquierda del cabello
  x      = rnorm(50)               
  y      = rnorm(50, mean = 2) 
  z      = t(cbind(t(x),t(y)));    # arma matriz (100x1)
  faces(as.matrix(z[1:50]),ind, main="Observaciones 1 a 50", ncol.plot = 5, nrow.plot = 10) # primeras 50 caras

## effect of variables:
##  modified item       Var   
##  "height of face   " "Var1"
##  "width of face    " "Var1"
##  "structure of face" "Var1"
##  "height of mouth  " "Var1"
##  "width of mouth   " "Var1"
##  "smiling          " "Var1"
##  "height of eyes   " "Var1"
##  "width of eyes    " "Var1"
##  "height of hair   " "Var1"
##  "width of hair   "  "Var1"
##  "style of hair   "  "Var1"
##  "height of nose  "  "Var1"
##  "width of nose   "  "Var1"
##  "width of ear    "  "Var1"
##  "height of ear   "  "Var1"
  faces(as.matrix(z[51:100]),ind, main="Observaciones 51 a 100", ncol.plot = 5, nrow.plot = 10) # primeras 50 caras

## effect of variables:
##  modified item       Var   
##  "height of face   " "Var1"
##  "width of face    " "Var1"
##  "structure of face" "Var1"
##  "height of mouth  " "Var1"
##  "width of mouth   " "Var1"
##  "smiling          " "Var1"
##  "height of eyes   " "Var1"
##  "width of eyes    " "Var1"
##  "height of hair   " "Var1"
##  "width of hair   "  "Var1"
##  "style of hair   "  "Var1"
##  "height of nose  "  "Var1"
##  "width of nose   "  "Var1"
##  "width of ear    "  "Var1"
##  "height of ear   "  "Var1"
  x
##  [1] -1.28507850  0.63355369  0.44602700  0.71122962  1.04101571  0.38166065
##  [7]  1.34567679  0.86589406 -0.35308261 -1.64286655 -0.20146989  1.27193157
## [13]  1.32968955 -0.52210557 -0.57437794  0.59510035  0.81025962 -0.46032635
## [19]  0.08005171  1.04435577  0.92007610  0.84796940 -1.25383572 -0.23593163
## [25] -0.52425121  0.69148768  0.66603702 -0.29236692  1.81754601  1.89572778
## [31] -0.02331735  0.16033679  0.16396680  0.30077960  0.34725328  0.15055643
## [37]  0.42060243 -0.05834590 -0.55080178 -1.23368553  0.21906266 -1.22663411
## [43] -0.47549157  1.11156992  0.02812435  1.92778317  0.23281437 -0.86684754
## [49] -2.15135470  0.84251688
  y
##  [1] 2.6937484 2.0669061 0.8156591 0.9635730 1.1769620 2.4175346 3.1611587
##  [8] 4.2057422 2.4187568 0.3011098 1.9511804 3.6623551 1.9038069 3.4798339
## [15] 1.4281953 0.3373708 1.5929995 0.3781713 1.7620376 2.6452088 0.8151953
## [22] 2.6671110 1.5816386 3.4534107 2.0263459 1.6859284 1.6761035 2.1650199
## [29] 2.0502183 3.0716692 1.5613290 1.9155633 2.8429410 1.3625479 2.4104710
## [36] 2.7806885 1.4568932 2.3856396 2.4258781 2.2510476 2.1942852 0.7077996
## [43] 3.5536185 3.3742319 2.1250455 2.9726510 1.4439805 1.6714803 0.8185543
## [50] 1.7004829

Sí hubo caras similares. Hay 4 payasos en Y comparado con los 3 que hay de X. Hay casi la misma cantidad de caras rojas y cabello verde; hay 4 caras amarillas en Y comparadas con las 10 que hay en X. Las caras rojas son aquellos valores cercanos a 0, mientras que los payasos son los valores más grandes (+/- 3)

3. Consideren los siguientes datos

a. Encuentra la proyección de \(X_1\) sobre 1’=(1,1,1,1,1,1)

  # x1 = Nómina de jugadores
  x1 <- c(3497900,2485475,1782875,1725450,1645575,1469800)

  # Defino al vector de 1 normalizado 
  ones <- rep(1,6)/sqrt(6)
  
  # Definición de proyección de x1 sobre 1
  proy = t(x1) %*% ones %*% ones
  proy
##         [,1]    [,2]    [,3]    [,4]    [,5]    [,6]
## [1,] 2101179 2101179 2101179 2101179 2101179 2101179

b. Calcula el vector desviación. Relaciona su longitud a la desviación estándar.

  # Definición de vector desviación
  vDesv <- x1 - proy
  vDesv
##         [,1]     [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 1396721 384295.8 -318304.2 -375729.2 -455604.2 -631379.2
  # Desviación estándar
  StanDev <- sd(x1)
  StanDev
## [1] 767752.2

Notemos que el cuadrado de la longitud del vector desviación es igual a la suma del cuadrado de las desviaciones

  # Longitud del vector
  norma <- norm(vDesv,type = "2")
  norma^2
## [1] 2.947217e+12
  vDesv%*%t(vDesv)
##              [,1]
## [1,] 2.947217e+12

c. Grafica (a escala) el triángulo formado por \(y_1\), \(\bar{x}_1\), \(y_1-\bar{x}_1\). Identifica la longitud de cada vector en tu gráfica

  nx1 = norm(x1, type = "2")
  nproy = norm(proy, type = "2")
  nvDesv = norm(vDesv, type = "2")

  
  x = c(0, nproy/nx1, nproy/nx1)
  y = c(0, 0,nvDesv/nx1)
  
  plot(x, y, xlim = c(0,1.1), ylim = c(0, .4))
  arrows(0,0, x1 = nproy/nx1, y1 = 0)
  arrows(0,0, x1 = nproy/nx1, y1 = nvDesv/nx1)
  arrows(nproy/nx1,0, x1 = nproy/nx1, y1 = nvDesv/nx1)
  
  vectores = cbind(x,y)
  text(vectores, labels = c(round(nproy/nx1,3),round(nvDesv/nx1,3),1), adj = c(0,-1))

d. Repetir los incisos (a) a (c) para \(X_2\) a’. Encuentra la proyección de \(X_2\) sobre 1’=(1,1,1,1,1,1)

  # x2 = % de perdidos-ganados
  x2 <- c(0.623,0.593,0.512,0.5,0.463,0.395)
  
  # Definición de proyección de x2 sobre 1
  proy2 = t(x2) %*% ones %*% ones
  proy2
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 0.5143333 0.5143333 0.5143333 0.5143333 0.5143333 0.5143333

b’. Calcula el vector de desviación

  # Definición de vector desviación
  vDesv2 <- x2 - proy2
  vDesv2
##           [,1]       [,2]         [,3]        [,4]        [,5]       [,6]
## [1,] 0.1086667 0.07866667 -0.002333333 -0.01433333 -0.05133333 -0.1193333
  # Desviación estándar
  StanDev2 <- sd(x2)
  StanDev2
## [1] 0.08376555

c’. Grafica (a escala) el triángulo formado por \(y_2\), \(\bar{x}_2\), \(y_2-\bar{x}_2\). Identifica la longitud de cada vector en tu gráfica

  nx2 = norm(x2, type = "2")
  nproy2 = norm(proy2, type = "2")
  nvDesv2 = norm(vDesv2, type = "2")

  
  x = c(0, nproy2/nx2, nproy2/nx2)
  y = c(0, 0,nvDesv2/nx2)
  
  plot(x, y, xlim = c(0,1.1), ylim = c(0, .4))
  arrows(0,0, x1 = nproy2/nx2, y1 = 0)
  arrows(0,0, x1 = nproy2/nx2, y1 = nvDesv2/nx2)
  arrows(nproy2/nx2,0, x1 = nproy2/nx2, y1 = nvDesv2/nx2)
  
  vectores = cbind(x,y)
  text(vectores, labels = c(round(nproy2/nx2,3),round(nvDesv2/nx2,3),1), adj = c(0,-1))

e. Grafica (a escala) los dos vectores desviación \(y_1-\bar{x}_1\) y \(y_2-\bar{x}_2\). Calcula el valor del ángulo entre ellos

  (theta<-acos(vDesv%*%t(vDesv2)/(nvDesv*nvDesv2)))
##           [,1]
## [1,] 0.4687649
  x = c(1,1*cos(theta))
  y = c(0,1*sin(theta))
  
  plot(x,y, xlim = c(0,2), ylim = c(0, 2))
  
  arrows(0,0, x1 = 1, y1 = 0)
  arrows(0,0,x1=1*cos(theta),y1=1*sin(theta) )

  vectores = cbind(x,y)

f. Calcula la varianza muestral generalizada (vmg) det(S) para estos datos e interpreta

  X <- cbind(x1,x2)                 # Creo la matriz de X1 y X2
  n <- nrow(X)
  uno <- rep(1,n)
  xBar <- t(X) %*% uno / n          # x barra
  H <- diag(n) - uno%*%t(uno) / n   # H
  Sn <- t(X) %*% H %*% X / n        # Sn
  S <- n/(n-1) * Sn
  
  vmg <- det(S)
  vmg
## [1] 844182191

Como la vmg es proporcional al cuadrado del olumen generado por los p vectores de desviación. COmo vmg = 844182191, entonces sabemos que el volumen del elipsoide al cuadrado generado por \(S\), será igual al vmg multiplicado por una constante proveniente de la ecuación del elipsoide correspondiente a \(S\).

g. Calcula la varianza muestral total (vmt) tr(S) para estos datos e interpreta

vmt <- sum(diag(S))
vmt
## [1] 589443426354

Geométricamente, la vmt es la suma de las longitudes al cuadrado de los p vectores de desviación divididos entre n-1; sin embargo, no le presta atención a la orientación de los vectores residuales. Como vmt = 589443426354, entonces sabemos que estas corresponde a la suma de las varrianzas de los p elementos de la matriz.

4. Dibuja las elipsoides sólidas para las tres matrices siguientes y determina los valores de los ejes mayores y menores

5. Archivo del INEGI

a. Configura la matriz X para poder operar con ella

# datos <- read.csv("D:\\ITAM\\Aplicada III\\Tareas\\HW2\\INEGIConstruccion2017.csv", header=FALSE)
datos <- read.csv(".\\INEGIConstruccion2017.csv", header=FALSE)
datos <- datos[-1,]

names(datos) <- as.matrix(datos[1,])    # Tomo primer renglón
datos <- datos[-1,]                     # Elimino el primer renglón
datos[] <- lapply(datos, function(x) type.convert(as.character(x))) # Primera columna names(datos) se vuelve el header

# Datos de enero 2017
datos <- datos %>% filter(Periodo == "2017/01")
head(datos)
# x1: total de horas trabajadas
horasTotalesEdos <- t(datos[,str_detect(names(datos), "Horas trabajadas > Total")])


#x2: valor total de producción
totalProd <- t(datos[,str_detect(names(datos), "tipo de obra > Total")])


#como faltan datos, falta agregar los "NA"s para que la matriz funcione correctamente.
#están los estados en orden y no falta ninguno antes de Michoacán, por lo que solamente hay que:
# 1) quitar el primer registro que es el total nacional
# 2) agregar 16 NAs hasta el final
totalProd <- as.matrix( totalProd[2:17,]) #1)
emptyVect <- as.matrix(rep(c(NA), times=16))
totalProd <- as.matrix(append(totalProd, emptyVect, after=16))


#x3: total de horas trabajadas Dependiente
horasDep <- as.matrix(t(datos[,str_detect(names(datos), "Dependiente > Total")])[-1,])


#x4: Obreros Dependiente
obrerosDep <- as.matrix(t(datos[,str_detect(names(datos), "Obreros")])[-1,])

#x5: Empleados Dependiente
empDep <- as.matrix(t(datos[,str_detect(names(datos), "Empleados")])[-1,])

#x6: Propietarios Dependiente
propDep <- as.matrix(t(datos[,str_detect(names(datos), "Propietarios")])[-1,])

#x7: total de horas trabajadas No dependiente
horasNoDep <- as.matrix(t(datos[,str_detect(names(datos), "No dependiente")])[-1,])

X <- cbind(horasTotalesEdos, totalProd, horasDep, obrerosDep, empDep, propDep, horasNoDep)

#Para hacer más clara nuestra tabla, adaptamos los nombres de los renglones
rownames(X) <- str_remove(as.array(rownames(X)), "Construcción \\(encuesta mensual\\) > Por entidad federativa > Horas trabajadas > Total ")
rownames(X) <- str_remove(as.array(rownames(X)), "\\(Miles de horas\\)")

X
##                                       [,1]      [,2]     [,3]     [,4]     [,5]
## Aguascalientes                    1441.083  725965.1 1353.470 1055.779  268.036
## Baja California                   3195.071  872017.8 2784.064 2045.011  723.914
## Baja California Sur                636.888  438205.6  564.851  329.325  200.202
## Campeche                          1669.002  533582.3 1323.013  958.880  349.655
## Coahuila de Zaragoza              6140.209 1059265.3 4001.635 3304.095  619.639
## Colima                            2189.691  414247.1 2152.042 1804.447  320.512
## Chiapas                           1469.275  717169.8 1422.475  893.406  492.414
## Chihuahua                         7164.411 1404806.6 6870.834 5722.790 1079.958
## Ciudad de México                 15689.435 2232549.7 9407.633 6836.521 2430.036
## Durango                           4352.566  599605.0 4303.292 3874.623  407.282
## Guanajuato                        5636.278 2559231.8 5493.818 4412.815  962.888
## Guerrero                          1511.229  345778.8 1435.678 1166.891  236.075
## Hidalgo                           1152.703  539917.0 1146.944  882.737  225.714
## Jalisco                           8089.720 2136522.1 7105.544 5830.045 1188.357
## México                            5392.139 3144478.7 5045.081 3118.509 1814.544
## Michoacán de Ocampo               4084.203  558029.3 3150.759 2363.378  715.655
## Morelos                            550.092        NA  511.352  356.827  129.513
## Nayarit                           1929.037        NA 1575.761 1273.005  271.445
## Nuevo León                       13389.724        NA 8528.908 6926.857 1530.889
## Oaxaca                             918.198        NA  881.406  549.736  249.908
## Puebla                            4487.444        NA 3833.660 3122.048  645.089
## Querétaro                         4037.281        NA 3778.835 3173.833  563.313
## Quintana Roo                      1947.345        NA 1765.056 1395.071  359.425
## San Luis Potosí                   2556.367        NA 2365.318 1799.165  511.367
## Sinaloa                           3345.227        NA 3143.382 2267.931  789.313
## Sonora                            6262.502        NA 5555.881 4631.798  911.309
## Tabasco                           2261.705        NA 1716.063 1169.822  501.094
## Tamaulipas                        6176.806        NA 5224.554 4208.895  950.857
## Tlaxcala                           268.718        NA  268.210  200.665   57.040
## Veracruz de Ignacio de la Llave   6582.393        NA 6226.744 5041.573 1113.837
## Yucatán                           3423.691        NA 3310.421 2572.381  708.740
## Zacatecas                         1510.317        NA 1495.929 1273.993  208.743
##                                     [,6]     [,7]
## Aguascalientes                    29.655   87.613
## Baja California                   15.139  411.007
## Baja California Sur               35.324   72.037
## Campeche                          14.478  345.989
## Coahuila de Zaragoza              77.901 2138.574
## Colima                            27.083   37.649
## Chiapas                           36.655   46.800
## Chihuahua                         68.086  293.577
## Ciudad de México                 141.076 6281.802
## Durango                           21.387   49.274
## Guanajuato                       118.115  142.460
## Guerrero                          32.712   75.551
## Hidalgo                           38.493    5.759
## Jalisco                           87.142  984.176
## México                           112.028  347.058
## Michoacán de Ocampo               71.726  933.444
## Morelos                           25.012   38.740
## Nayarit                           31.311  353.276
## Nuevo León                        71.162 4860.816
## Oaxaca                            81.762   36.792
## Puebla                            66.523  653.784
## Querétaro                         41.689  258.446
## Quintana Roo                      10.560  182.289
## San Luis Potosí                   54.786  191.049
## Sinaloa                           86.138  201.845
## Sonora                            12.774  706.621
## Tabasco                           45.147  545.642
## Tamaulipas                        64.802  952.252
## Tlaxcala                          10.505    0.508
## Veracruz de Ignacio de la Llave   71.334  355.649
## Yucatán                           29.300  113.270
## Zacatecas                         13.193   14.388

b. Calcula el vector de medias y la matriz de covarianzas de la matrix X

(XBarra <- colMeans(X, na.rm = TRUE))
## [1] 4.045648e+03 1.142586e+06 3.366957e+03 2.642589e+03 6.730238e+02
## [6] 5.134369e+01 6.786918e+02
(S <- var(X, na.rm = TRUE))
##              [,1]         [,2]         [,3]         [,4]         [,5]
## [1,]   14567990.6   2232551610 9.412867e+06 7.148303e+06   2141883.71
## [2,] 2232551610.2 781734140931 1.730687e+09 1.234355e+09 465994225.21
## [3,]    9412867.2   1730687190 6.726512e+06 5.223303e+06   1420140.80
## [4,]    7148303.3   1234355219 5.223303e+06 4.159499e+06   1003797.58
## [5,]    2141883.7    465994225 1.420141e+06 1.003798e+06    394867.16
## [6,]     122680.2     30337746 8.306814e+04 6.000603e+04     21476.06
## [7,]    5155123.4    501864420 2.686356e+06 1.925001e+06    721742.90
##              [,6]         [,7]
## [1,]   122680.183   5155123.42
## [2,] 30337745.581 501864419.99
## [3,]    83068.142   2686355.70
## [4,]    60006.032   1925000.76
## [5,]    21476.064    721742.90
## [6,]     1586.046     39612.04
## [7,]    39612.041   2468767.72

c. Verifiquen que det(S) = \(\prod_{i = 1}^{7} s_{i}^2\) \(\times\)det(R)

Si analizamos los eigenvalores de nuestra matriz de varianza de los datos, notamos que los valores son muy dispersos, y uno de ellos es muy cercano a cero. Esto puede ocasionar que su producto sea extremádamente cercano a cero y el resultado deseado no se alcance. Por ende, vamos a aplicar la función \(log()\) a los datos, para corregir así la disparidad de los eigenvalores.

X <- log(X)
R <- cor(X, use = "complete.obs")
detR <- det(R)
S <- var(X, na.rm = TRUE)
detS <- det(S)

vaps <- eigen(S)$values

(det(prod(vaps)*cor(R))-det(S))
## [1] -9.496053e-09

Como podemos notar, la diferencia es casi cero, por lo que podemos decir que el error es despreciable y así afirmamos que el enunciado de arriba es cierto.

Pregunta 6

Pregunta 7

Pregunta 8

Pregunta 9